!############################# Change Log ##################################
! 5.0.0
!
!###########################################################################
!  Copyright (C)  1990, 1995, 1999, 2000, 2003 - All Rights Reserved
!  Regional Atmospheric Modeling System - RAMS
!###########################################################################

subroutine chem_makevarf(ng)

use isan_coms

!srf-chem
use chem_isan_coms
use mem_chem1, only:  CHEM_ASSIM     ! intent(in)
!srf-chem-end

use mem_grid

implicit none

integer :: ng

!---------------------------------------------------------------+
!    Interpolate model grid from isentropic/sigmaz/surface data
!---------------------------------------------------------------+

!            Vertically interpolate isentropic data to sigma-z levels

call isnsig(nnzp(ng),nnxp(ng),nnyp(ng) ,is_grids(ng)%rr_u(1,1,1)  &
        ,is_grids(ng)%rr_v(1,1,1)      ,is_grids(ng)%rr_t(1,1,1)  &
        ,is_grids(ng)%rr_r(1,1,1)      ,is_grids(ng)%rr_p(1,1,1)  &
        ,grid_g(ng)%topt(1,1),ztn(1,ng),ztop)

!srf-chem
if(CHEM_ASSIM == 1 .and. nspecies>0 ) then
    call chem_isnsig(nnzp(ng),nnxp(ng),nnyp(ng),nspecies     &
                    ,chem_is_grids(ng)%rr_sc(1,1,1,1)        &
                    ,grid_g(ng)%topt(1,1),ztn(1,ng),ztop)
endif
!srf-chem-end

!            Compute Exner function on model sigma-z surfaces
!              and change relative humidity to mixing ratio.

call vshyd(nnzp(ng),nnxp(ng),nnyp(ng),is_grids(ng)%rr_p(1,1,1)  &
     ,is_grids(ng)%rr_t(1,1,1)       ,is_grids(ng)%rr_r(1,1,1)  &
     ,grid_g(ng)%topt(1,1),grid_g(ng)%rtgt(1,1),ztn(1,ng))

!          Combine surface analysis with the upper air data.

call visurf(nnzp(ng),nnxp(ng),nnyp(ng) ,is_grids(ng)%rr_u(1,1,1)  &
      ,is_grids(ng)%rr_v(1,1,1)        ,is_grids(ng)%rr_t(1,1,1)  &
      ,is_grids(ng)%rr_r(1,1,1)        ,is_grids(ng)%rr_p(1,1,1)  &
      ,grid_g(ng)%topt(1,1),grid_g(ng)%rtgt(1,1),ztn(1,ng))

is_grids(ng)%rr_slp (1:nnxp(ng),1:nnyp(ng)) =rs_slp (1:nnxp(ng),1:nnyp(ng))
is_grids(ng)%rr_sfp (1:nnxp(ng),1:nnyp(ng)) =rs_sfp (1:nnxp(ng),1:nnyp(ng))
is_grids(ng)%rr_sft (1:nnxp(ng),1:nnyp(ng)) =rs_sft (1:nnxp(ng),1:nnyp(ng))
is_grids(ng)%rr_snow(1:nnxp(ng),1:nnyp(ng)) =rs_snow(1:nnxp(ng),1:nnyp(ng))
is_grids(ng)%rr_sst (1:nnxp(ng),1:nnyp(ng)) =rs_sst (1:nnxp(ng),1:nnyp(ng))

!          average the velocities to the correct points in the stagger
!             and rotate for polar stereographic transformation.

call varuv(nnzp(ng),nnxp(ng),nnyp(ng),is_grids(ng)%rr_u(1,1,1)  &
                      ,is_grids(ng)%rr_v(1,1,1))

return
end
!
!     **************************************************************
!
!subroutine isnsig(n1,n2,n3,uu,vv,tt,rr,pp  &
!     ,topt,zt,ztop)
!   
!use isan_coms
!use rconstants
!
!implicit none
!
!integer :: n1,n2,n3
!real, dimension(n1,n2,n3) :: uu,vv,tt,rr,pp
!real, dimension(n2,n3) :: topt
!real, dimension(n1) :: zt
!real :: ztop
!
!real, allocatable :: v1(:),v2(:),v3(:),v4(:),v5(:),v6(:)
!integer :: i,j,k,nki,ki,kbeg,kend
!real :: rtg,wtsz
!
!allocate (v1(n1),v2(nisn),v3(nisn),v4(nisn),v5(nisn),v6(nisn))
!
!do j=1,n3
!   do i=1,n2
!
!      rtg=(1.-topt(i,j)/ztop)
!      do k=1,n1
!	  v1(k)=topt(i,j)+zt(k)*rtg
!      enddo
!
!      if (guess1st.eq.'PRESS') then
!	  nki=0
!	  do ki=1,nisn
!!	       print*,'v0:',i,j,ki,pi_s(i,j,ki),levth(ki),pi_p(i,j,ki)  &
!!		 ,pi_u(i,j,ki),pi_v(i,j,ki)
!	     if(pi_p(i,j,ki).lt.1e20.and.pi_s(i,j,ki).lt.1e20  &
!		  .and.pi_u(i,j,ki).lt.1e20.and.pi_v(i,j,ki).lt.1e20)  &
!		  then
!		nki=nki+1
!		v2(nki)=(pi_s(i,j,ki)-cp*levth(ki)  &
!		     *(pi_p(i,j,ki)*p00i)**rocp)/g
!!		      print*,'v2:',i,j,nki,pi_s(i,j,ki),levth(ki),pi_p(i,j,ki)
!		     
!		v3(nki)=pi_u(i,j,ki)
!		v4(nki)=pi_v(i,j,ki)
!		v5(nki)=levth(ki)
!		v6(nki)=pi_r(i,j,ki)
!	     endif
!	  enddo
!
!	  call htint(nki,v3,v2,n1,uu(1,i,j),v1)
!	  call htint(nki,v4,v2,n1,vv(1,i,j),v1)
!	  call htint(nki,v5,v2,n1,tt(1,i,j),v1)
!	  call htint(nki,v6,v2,n1,rr(1,i,j),v1)
!
!      endif
!
!
!      kbeg=n1+1
!      kend=n1+1
!      do k=1,nsigz
!	  if(zt(k).gt.hybbot) then
!	     kbeg=k-1
!	     exit
!	  endif
!      enddo
!
!      do k=1,nsigz
!	  if(zt(k).gt.hybtop) then
!	     kend=k
!	     exit
!	  endif
!      enddo
!
!      do k=1,nsigz
!
!	  if(k.lt.kbeg) then
!	     wtsz=sigzwt
!	  elseif(k.gt.kend) then
!	     wtsz=0.
!	  elseif(k.ge.kbeg.and.k.le.kend) then
!	     wtsz= (zt(kend)-zt(k))  &
!		  /(zt(kend)-zt(kbeg))  &
!		  *sigzwt
!	  endif
!
!	  uu(k,i,j)=(1.-wtsz)*uu(k,i,j)+wtsz*ps_u(i,j,k)
!	  vv(k,i,j)=(1.-wtsz)*vv(k,i,j)+wtsz*ps_v(i,j,k)
!	  tt(k,i,j)=(1.-wtsz)*tt(k,i,j)+wtsz*ps_t(i,j,k)
!	  rr(k,i,j)=(1.-wtsz)*rr(k,i,j)+wtsz*ps_r(i,j,k)
!
!	  if(guess1st.eq.'RAMS') then
!	     pp(k,i,j)=ps_p(i,j,k)
!	  endif
!
!      !   if(i==10.and.j==10) &
!      !      print*,'avrf, k:',k,tt(k,i,j),rr(k,i,j), pp(k,i,j)
!
!      enddo
!
!   enddo
!enddo
!
!deallocate (v1,v2,v3,v4,v5,v6)
!
!return
!end
!
!!     ****************************************************************
!
!subroutine visurf(n1,n2,n3,up,vp,thp,rtp,pp,topt,rtgt,zt)
!
!use isan_coms
!use rconstants
!
!implicit none
!
!integer :: n1,n2,n3
!real, dimension(n1,n2,n3) :: up,vp,pp,thp,rtp
!real, dimension(n2,n3) :: rtgt,topt
!real, dimension(n1) :: zt
!
!integer :: i,j,k 
!
!real, allocatable :: v3(:)
!real :: zgp,dzgobs,wt,ppp,ttt,rhm
!real, external :: rs
!
!
!allocate(v3(n1))
!
!do j=1,n3
!   do i=1,n2
!      if(rs_qual(i,j).gt..5) then
!
!!	    Quality of this surface analysis point is okay.
!
!!	    Extend the vertical influence of the surface analysis
!!	      up to a distance of SFCINF above the height of the
!!	      objectively analyzed surface.
!
!	  do k=1,n1
!	     zgp=topt(i,j)+zt(k)*rtgt(i,j)
!	     dzgobs=abs(zgp-rs_top(i,j))
!	     wt=dzgobs/sfcinf
!	     if(wt.le.1.) then
!		if(rs_u(i,j).lt.1.e10)  &
!		     up(k,i,j)=up(k,i,j)*wt+rs_u(i,j)*(1.-wt)
!		if(rs_v(i,j).lt.1.e10)  &
!		     vp(k,i,j)=vp(k,i,j)*wt+rs_v(i,j)*(1.-wt)
!		ppp=(pp(k,i,j)/cp)**cpor*p00
!		ttt=thp(k,i,j)*pp(k,i,j)/cp
!		rhm=rtp(k,i,j)/rs(ppp,ttt)
!		if(rs_r(i,j).lt.1.e10)  &
!		     rhm=rhm*wt+rs_r(i,j)*(1.-wt)
!		if(rs_t(i,j).lt.1.e10)  &
!		     thp(k,i,j)=thp(k,i,j)*wt+rs_t(i,j)*(1.-wt)
!		ttt=thp(k,i,j)*pp(k,i,j)/cp
!		rtp(k,i,j)=rhm*rs(ppp,ttt)
!	     endif
!	  enddo
!
!	  do k=1,n1
!	     v3(k)=topt(i,j)+zt(k)*rtgt(i,j)
!	  enddo
!	  do k=n1-1,1,-1
!	     pp(k,i,j)=pp(k+1,i,j)+g*(v3(k+1)-v3(k))  &
!		  /((thp(k,i,j)*(1.+.61*rtp(k,i,j))+thp(k+1,i,j)  &
!		  *(1.+.61*rtp(k+1,i,j)))*.5)
!	  enddo
!
!      endif
!
!   enddo
!enddo
!
!deallocate(v3)
!
!return
!end
!
!!     ****************************************************************
!
!subroutine vshyd(n1,n2,n3,pp,tt,rr,topt,rtg,zt)
!     
!use isan_coms
!use rconstants
!
!implicit none
!     
!integer :: n1,n2,n3
!real, dimension(n1,n2,n3) :: tt,rr,pp
!real, dimension(n2,n3) :: topt,rtg
!real, dimension(n1) :: zt
!
!integer :: i,j,k,lbc,kabc
!real :: thmin,ppp,ttt,tmpbc,rvibc,piibc,ziibc,gd2
!real, external :: rs
!
!real, allocatable :: v3(:)
!
!allocate(v3(n1))
!
!if(guess1st.eq.'PRESS') then
!
!!	   Find internal boundary condition for the hydrostatic
!!	     equation
!
!   thmin=tt(n1-1,1,1)
!   do j=1,n3
!      do i=1,n2
!	  thmin=min(thmin,tt(n1-1,i,j))
!      enddo
!   enddo
!
!   do k=nisn,1,-1
!      print*,'thmins-',k,levth(k),thmin
!      if(float(levth(k)).le.thmin) then
!	  lbc=k
!	  go to 121
!      endif
!   enddo
!   stop 'vi-nobc'
!121	 continue
!
!   print 1212,lbc,levth(lbc)
!1212	 format(//' Hydrostatic boundary condition set at level'  &
!	 ,2I4,' K')
!
!endif
!
!
!
!do j=1,n3
!   do i=1,n2
!
!!	   Compute actual model heights
!
!      do k=1,n1
!	  v3(k)=topt(i,j)+zt(k)*rtg(i,j)
!      enddo
!
!      if(guess1st.eq.'PRESS') then
!
!!	   Find model level above boundary condition
!
!	  do k=1,n1
!	     if(tt(k,i,j).gt.levth(lbc))go to 141
!	  enddo
!	  stop 'vi-nbc2'
!141	  continue
!	  kabc=k
!
!!	   Integrate P to surface
!
!	  piibc=cp*(pi_p(i,j,lbc)*p00i)**rocp
!	  ziibc=(pi_s(i,j,lbc)-levth(lbc)*piibc)/g
!	  gd2=2.*g
!	  pp(kabc-1,i,j)=piibc+gd2*(ziibc-v3(kabc-1))  &
!	       /(levth(lbc)+tt(kabc-1,i,j))
!	  do k=kabc-2,1,-1
!	     pp(k,i,j)=pp(k+1,i,j)+gd2*(v3(k+1)-v3(k))  &
!		  /(tt(k+1,i,j)+tt(k,i,j))
!	  enddo
!
!!	   Integrate P to model top
!
!	  do k=kabc,n1
!	     pp(k,i,j)=pp(k-1,i,j)-gd2*(v3(k)-v3(k-1))  &
!		  /(tt(k,i,j)+tt(k-1,i,j))
!	  enddo
!!
!!	   Compute mixing ratio from relative humidity and do final
!!	     hydrostatic integration.
!!
!	  do k=1,n1
!	     ppp=(pp(k,i,j)/cp)**cpor*p00
!	     ttt=tt(k,i,j)*pp(k,i,j)/cp
!	     rr(k,i,j)=rr(k,i,j)*rs(ppp,ttt)
!	  enddo
!
!	  tmpbc=levth(lbc)*piibc/cp
!	  rvibc=pi_r(i,j,lbc)*rs(pi_p(i,j,lbc),tmpbc)
!
!	  pp(kabc-1,i,j)=piibc+gd2*(ziibc-v3(kabc-1))  &
!	       /(levth(lbc)*(1.+.61*rvibc)  &
!		+tt(kabc-1,i,j)*(1.+.61*rr(kabc-1,i,j)))
!	  do k=kabc-2,1,-1
!	     pp(k,i,j)=pp(k+1,i,j)+gd2*(v3(k+1)-v3(k))  &
!		  /(tt(k+1,i,j)*(1.+.61*rr(k+1,i,j))  &
!		   +tt(k,i,j)*(1.+.61*rr(k,i,j)))
!	  enddo
!
!	  do k=kabc,n1
!	     pp(k,i,j)=pp(k-1,i,j)-gd2*(v3(k)-v3(k-1))  &
!		  /((tt(k,i,j)*(1.+.61*rr(k,i,j))+tt(k-1,i,j)  &
!		  *(1.+.61*rr(k-1,i,j))))
!	  enddo
!
!      elseif (guess1st.eq.'RAMS') then
!
!	  do k=1,n1
!	     ppp=pp(k,i,j)
!	     ttt=tt(k,i,j)*(ppp/p00)**rocp
!	     rr(k,i,j)=rr(k,i,j)*rs(ppp,ttt)
!	  enddo
!
!	  pp(1,i,j)= cp*(pp(1,i,j)/p00)**rocp
!	  do k=2,n1
!	     pp(k,i,j)=pp(k-1,i,j)-g*(v3(k)-v3(k-1))  &
!		  /((tt(k,i,j)*(1.+.61*rr(k,i,j))+tt(k-1,i,j)  &
!		  *(1.+.61*rr(k-1,i,j)))*.5)
!	!  if(i==10.and.j==10) &
!	!     print*,'avrf222, k:',k,tt(k,i,j),rr(k,i,j), pp(k,i,j)
!	  enddo
!
!      endif
!
!   enddo
!enddo
!
!deallocate(v3)
!
!return
!end
!
!!     ******************************************************************
!
!subroutine varuv(n1,n2,n3,up,vp)
!
!implicit none
!
!integer :: n1,n2,n3
!real, dimension(n1,n2,n3) :: up,vp
!
!integer :: i,j,k
!
!!	       Average to u,v points
!
!do j=1,n3
!   do i=1,n2-1
!      do k=1,n1
!	  up(k,i,j)=(up(k,i,j)+up(k,i+1,j))*.5
!      enddo
!   enddo
!enddo
!
!do j=1,n3-1
!   do i=1,n2
!      do k=1,n1
!	  vp(k,i,j)=(vp(k,i,j)+vp(k,i,j+1))*.5
!      enddo
!   enddo
!enddo
!
!return
!end
!!
!
!     ****************************************************************

subroutine chem_varfile_nstfeed(ifm,icm,n1f,n2f,n3f,n1c,n2c,n3c &
                               ,nbot,ntop)

use isan_coms
!srf-chem
use chem_isan_coms
use mem_chem1, only:  CHEM_ASSIM     ! intent(in)
!srf-chem-end

implicit none

integer :: ifm,icm,n1f,n2f,n3f,n1c,n2c,n3c,nbot,ntop,nspc


!     Feed back the finer mesh to the coarser mesh.

call fdback(is_grids(icm)%rr_u   (1,1,1),is_grids(ifm)%rr_u   (1,1,1) &
           ,is_grids(icm)%rr_dn0u(1,1,1),is_grids(ifm)%rr_dn0u(1,1,1) &
           ,n1c,n2c,n3c,n1f,n2f,n3f,ifm,'u',rr_scr1(1))

call fdback(is_grids(icm)%rr_v   (1,1,1),is_grids(ifm)%rr_v   (1,1,1) &
           ,is_grids(icm)%rr_dn0v(1,1,1),is_grids(ifm)%rr_dn0v(1,1,1) &
           ,n1c,n2c,n3c,n1f,n2f,n3f,ifm,'v',rr_scr1(1))

call fdback(is_grids(icm)%rr_p  (1,1,1),is_grids(ifm)%rr_p  (1,1,1) &
           ,is_grids(icm)%rr_dn0(1,1,1),is_grids(ifm)%rr_dn0(1,1,1) &
           ,n1c,n2c,n3c,n1f,n2f,n3f,ifm,'p',rr_scr1(1))

call fdback(is_grids(icm)%rr_t  (1,1,1),is_grids(ifm)%rr_t  (1,1,1) &
           ,is_grids(icm)%rr_dn0(1,1,1),is_grids(ifm)%rr_dn0(1,1,1) &
           ,n1c,n2c,n3c,n1f,n2f,n3f,ifm,'t',rr_scr1(1))

call fdback(is_grids(icm)%rr_r  (1,1,1),is_grids(ifm)%rr_r  (1,1,1) &
           ,is_grids(icm)%rr_dn0(1,1,1),is_grids(ifm)%rr_dn0(1,1,1) &
           ,n1c,n2c,n3c,n1f,n2f,n3f,ifm,'t',rr_scr1(1))



!srf-chem
if(CHEM_ASSIM == 1 .and. nspecies>0 ) then
 do nspc=1,nspecies
  print*,'fdback for spc=',nspc
  call fdback(chem_is_grids(icm)%rr_sc  (1,1,1,nspc) &
             ,chem_is_grids(ifm)%rr_sc  (1,1,1,nspc) &
             ,is_grids(icm)%rr_dn0(1,1,1),is_grids(ifm)%rr_dn0(1,1,1) &
             ,n1c,n2c,n3c,n1f,n2f,n3f,ifm,'t',rr_scr1(1))
 enddo
endif
!srf-chem-end



if(nbot == 1) then
   call botset(n1c,n2c,n3c,1,n2c,1,n3c  &
        ,15 ,is_grids(icm)%rr_u(1,1,1),'U')
   call botset(n1c,n2c,n3c,1,n2c,1,n3c  &
        ,15 ,is_grids(icm)%rr_v(1,1,1),'V')
   call botset(n1c,n2c,n3c,1,n2c,1,n3c  &
        ,15 ,is_grids(icm)%rr_p(1,1,1),'P')
   call botset(n1c,n2c,n3c,1,n2c,1,n3c  &
        ,15 ,is_grids(icm)%rr_t(1,1,1),'T')
   call botset(n1c,n2c,n3c,1,n2c,1,n3c  &
        ,15 ,is_grids(icm)%rr_r(1,1,1),'T')
endif

if(ntop == 1) then
   call topset(n1c,n2c,n3c,1,n2c,1,n3c  &
        ,15,is_grids(icm)%rr_u(1,1,1),is_grids(icm)%rr_u(1,1,1),'U')
   call topset(n1c,n2c,n3c,1,n2c,1,n3c  &
        ,15,is_grids(icm)%rr_v(1,1,1),is_grids(icm)%rr_v(1,1,1),'V')
   call topset(n1c,n2c,n3c,1,n2c,1,n3c  &
        ,15,is_grids(icm)%rr_p(1,1,1),is_grids(icm)%rr_p(1,1,1),'P')
   call topset(n1c,n2c,n3c,1,n2c,1,n3c  &
        ,15,is_grids(icm)%rr_t(1,1,1),is_grids(icm)%rr_t(1,1,1),'T')
   call topset(n1c,n2c,n3c,1,n2c,1,n3c  &
        ,15,is_grids(icm)%rr_r(1,1,1),is_grids(icm)%rr_r(1,1,1),'T')
endif


return
end
!     **************************************************************
!
!srf-chem
subroutine chem_isnsig(n1,n2,n3,nsp_dummy,sc,topt,zt,ztop)
   
use isan_coms

!srf-chem
use chem_isan_coms
!srf-chem-end

use rconstants

implicit none

!srf-chem
integer :: n1,n2,n3,nsp_dummy
real, dimension(n1,n2,n3,nsp_dummy) :: sc
!srf-chem-end
real, dimension(n2,n3) :: topt
real, dimension(n1) :: zt
real :: ztop

!srf-chem
real, allocatable :: v1(:),v2(:),v6(:,:)
integer :: i,j,k,nki,ki,kbeg,kend,nspc
!srf-chem-end
real :: rtg,wtsz

!srf-chem
allocate (v1(n1),v2(nisn),v6(nisn,nsp_dummy))
!srf-chem-end

do j=1,n3
   do i=1,n2

      rtg=(1.-topt(i,j)/ztop)
      do k=1,n1
         v1(k)=topt(i,j)+zt(k)*rtg
      enddo

      if (guess1st.eq.'PRESS') then
!srf-chem
         do nspc=1,nsp_dummy
            nki=0
            do ki=1,nisn
              if(pi_p(i,j,ki).lt.1e20 .and. pi_s(i,j,ki).lt.1e20  &
                 .and. pi_sc(i,j,ki,nspc).lt.1e20)  then
                 nki=nki+1
                 v2(nki)=(pi_s(i,j,ki)-cp*levth(ki)  &
                      *(pi_p(i,j,ki)*p00i)**rocp)/g
                 v6(nki,nspc)=pi_sc(i,j,ki,nspc)
              endif
            enddo  ! loop ki
         enddo   ! loop nspc

!srf-chem
         do nspc=1,nsp_dummy
          call htint(nki,v6(1,nspc),v2,n1,sc(1,i,j,nspc),v1)
         enddo  ! loop nspc

      endif


      kbeg=n1+1
      kend=n1+1
      do k=1,nsigz
         if(zt(k).gt.hybbot) then
            kbeg=k-1
            exit
         endif
      enddo

      do k=1,nsigz
         if(zt(k).gt.hybtop) then
            kend=k
            exit
         endif
      enddo

      do k=1,nsigz

         if(k.lt.kbeg) then
            wtsz=sigzwt
         elseif(k.gt.kend) then
            wtsz=0.
         elseif(k.ge.kbeg.and.k.le.kend) then
            wtsz= (zt(kend)-zt(k))  &
                 /(zt(kend)-zt(kbeg))  &
                 *sigzwt
         endif

!srf-chem
         do nspc=1,nsp_dummy
           sc(k,i,j,nspc)=(1.-wtsz)*sc(k,i,j,nspc)+wtsz*ps_sc(i,j,k,nspc)
         enddo  ! loop nspc

      enddo

   enddo
enddo

deallocate (v1,v2,v6)

return
end subroutine chem_isnsig
!srf-chem-end
!     ****************************************************************

